! This is a 1D advection example using square initial condition and periodic
! boundary condition for MPDATA finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-21: Initial creation.
! - 2024-03-24: Use adv_1d_test_case_mod.

program mpdata_adv_1d_case

  use adv_1d_test_case_mod

  implicit none

  real, allocatable :: rho(:,:)     ! Tracer density being advected at cell centers
  real, allocatable :: flx(:)       ! Flux at cell interfaces
  real, allocatable :: uc (:)       ! Antidiffusion velocity
  integer :: iord = 3               ! Scheme order
  real, parameter :: eps = 1.0d-15  ! A small value to avoid divided-by-zero
  integer, parameter :: ns = 1      ! Stencil width
  integer i
  character(30), parameter :: scheme = 'mpdata'

  namelist /params/ nx, nt, dt, iord, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  allocate(flx(1-ns:nx+ns))
  allocate(uc (1-ns:nx+ns))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call mpdata(rho(:,old), rho(:,new), flx)
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  deallocate(rho, flx, uc)

  call adv_1d_test_case_final()

contains

  subroutine mpdata(q, qstar, flx)

    real, intent(in   ) :: q    (1-ns:nx+ns)
    real, intent(inout) :: qstar(1-ns:nx+ns)
    real, intent(out  ) :: flx  (1-ns:nx+ns)

    real c
    integer i, j

    c = dt / dx
    uc(:) = u
    qstar(:) = q(:)
    do j = 1, iord
      if (j > 1) then
        ! Calculate antidiffusion velocity.
        do i = 1, nx
          uc(i) = (abs(uc(i)) - uc(i)**2 * c) * &
                  (qstar(i+1) - qstar(i)) / (qstar(i+1) + qstar(i) + eps)
        end do
      end if
      call upwind(qstar, uc, flx)
      ! Update tracer density.
      do i = 1, nx
        qstar(i) = qstar(i) - c * (flx(i) - flx(i-1))
      end do
      call apply_bc(ns, nx, qstar)
    end do

  end subroutine mpdata

  subroutine upwind(q, u, f)

    real, intent(in ) :: q(1-ns:nx+ns)
    real, intent(in ) :: u(1-ns:nx+ns)
    real, intent(out) :: f(1-ns:nx+ns)

    integer i

    do i = 1, nx
      flx(i) = 0.5d0 * (u(i) * (q(i+1) + q(i)) - abs(u(i)) * (q(i+1) - q(i)))
    end do
    call apply_bc(ns, nx, flx)

  end subroutine upwind

end program mpdata_adv_1d_case
